disp('fmdlf');

h0 = 0.05;




fprintf('max loss = %e \n', max(abs(v_diss)));

% close(figure(6));

disp('fmdlf');

Num = 200;
vy2 = fn_total_solve(Num, h0);
v_diss2 = fn_v_diss(vy2, h0);

Num = 400;
vy4 = fn_total_solve(Num, h0);
v_diss4 = fn_v_diss(vy4, h0);

Num = 600;
vy6 = fn_total_solve(Num, h0);
v_diss6 = fn_v_diss(vy6, h0);

Num = 800;
vy8 = fn_total_solve(Num, h0);
v_diss8 = fn_v_diss(vy8, h0);



fprintf('the error of %d : %e \n', 2, max(abs(v_diss2)));
fprintf('the error of %d : %e \n', 4, max(abs(v_diss4)));
fprintf('the error of %d : %e \n', 6, max(abs(v_diss6)));
fprintf('the error of %d : %e \n', 8, max(abs(v_diss8)));

close(figure(99));
figure(99);
hold on;
plot(1:length(vy2), vy2, 'o');
plot(1:length(vy4), vy4, 'x');
plot(1:length(vy6), vy6, 's');
plot(1:length(vy8), vy8, '*');

% plot(1:length(vy8), h0*(1:length(vy8)), '-bo');



disp('fmdkx');

% Num = 1000;
% h0 = 0.;
% 
% [s_left, s_right] = fn_total_solve_lr(Num, h0);
% vy = fn_total_solve(Num, h0);
% 
% v_diss = zeros(1, length(vy));
% 
% 
% for n = 1:length(vy)
%     
%    tmp = vy(n) + f(vy(n));
%    for n1 = 1:(n-1)
%        tmp = tmp + 2*f(vy(n1));
%    end
%    
%    tmp = tmp - n*h0 - 2*n*sum(log(1+exp(-vy)));
% %    fprintf('n = %d, diss = %e \n', n, tmp);
%    v_diss(n) = tmp;
%     
% end
% 
% fprintf('max diss = %e \n', max(abs(v_diss)));
% figure(1); hold on;
% plot(1:length(vy), vy, '-o');

vh = 0.1:0.1:2.8;
vb = zeros(1, length(vh));
vsum = zeros(1, length(vh));
for ih = 1:length(vh)
    h0 = vh(ih);
    vy = fn_total_solve(Num, h0);
    vb(ih) = fn_b(vy, h0);
    vsum(ih) = sum(log(1+exp(-vy)));
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function b = fn_b(vy,h0)

b = vy(end) - length(vy)*h0;

end


function v_diss = fn_v_diss(vy, h0)

v_diss = zeros(1, length(vy));

f = @(x) log(1+exp(-x));
for n = 1:length(vy)
    
   tmp = vy(n) + f(vy(n));
   for n1 = 1:(n-1)
       tmp = tmp + 2*(n - n1)*f(vy(n1));
   end
   
   tmp = tmp - n*h0 - 2*n*sum(log(1+exp(-vy)));
%    fprintf('n = %d, diss = %e \n', n, tmp);
   tmp_real = real(tmp);
   tmp_imag = imag(tmp);
   tmp_imag = mod(tmp_imag+pi, 2*pi) - pi;
   tmp = tmp_real + 1j*tmp_imag;
   v_diss(n) = tmp;
    
end

end


function res = fn_find_singularity(Num, h0)

s_left = 0.1;
while ~fn_if_singu(Num, h0, s_left)
    s_left = s_left - 0.1;
end

s_right = 0.9;
while fn_if_singu(Num, h0, s_right)
    s_right = s_right + 0.1;
end

while 1
    
    if s_right - s_left < 1e-12
        break;
    end
    
    s_mid = 1/2*(s_left + s_right);
    if fn_if_singu(Num, h0, s_mid)
        s_left = s_mid;
    else
        s_right = s_mid;
    end

end

res = s_right;

end




function if_singu = fn_if_singu(Num, h0, s0)

vy = fn_solve_vy(Num, h0, s0);
tmp = sum(abs(imag(vy)));
if tmp>0.1
    if_singu = (1>0);
else
    if_singu = (1<0);
end

end



function [vy, s_mid] = fn_total_solve(Num, h0)

[s_left, s_right] = fn_total_solve_lr(Num, h0);

if fn_judge(Num, h0, s_left) == 0
    vy = fn_solve_vy(Num, h0, s_left);
    return;
end

if fn_judge(Num, h0, s_right) == 0
    vy = fn_solve_vy(Num, h0, s_right);
    return;
end

while 1

    s_mid = 1/2*(s_left + s_right);
    
    if abs(s_left - s_right) < 1e-14
        break;
    end
    
    if fn_judge(Num, h0, s_mid) == 0
        vy = fn_solve_vy(Num, h0, s_mid);
        return;
    elseif fn_judge(Num, h0, s_mid)>0
        s_right = s_mid;
    else
        s_left = s_mid;
    end
    
end

vy = fn_solve_vy(Num, h0, s_mid);

end



function [s_left, s_right] = fn_total_solve_lr(Num, h0)

s_left = fn_find_singularity(Num, h0);

s_right = s_left;

while 1
    s_right = s_right + 0.1;
    if fn_judge(Num, h0, s_right) > 0
        break;
    end
end

s_left = s_right - 0.1;

fprintf('s_left = %f, s_right = %f \n', s_left, s_right);
fprintf('judge_left = %f, judge_right = %f, \n', ...
    fn_judge(Num, h0, s_left), fn_judge(Num, h0, s_right));

end



function res = fn_judge(Num, h0, s)

res = s - real(fn_s_to_s(Num, h0, s));

end


function exps = fn_s_to_exps(Num, h0, s0)

vy = fn_solve_vy(Num, h0, s0);
exps = fn_exps(vy);

end


function [s, vy] = fn_s_to_s(Num, h0, s0)

vy = fn_solve_vy(Num, h0, s0);
s = fn_s(vy);

end


function exps = fn_exps(vy)

exps = prod(1 + exp(-vy));

exps = real(exps);

end


function s = fn_s(vy)

s = sum(log(1 + exp(-vy)));

end


function vy = fn_solve_vy(Num, h0, s)

W = h0 + 2*s;
vy = zeros(1, Num);

for n = 1:Num
    
    tv = vy(1:n-1);
    RHS = n*W - 2*sum(((n-1):-1:1).*log(1 + exp(-tv)));
    vy(n) = fn_solve(RHS);
    
end



end



function res = fn_solve(f)

if f < 33
    res = log(exp(f) - 1);
else
    res = f;
end


end